home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / complex.1 < prev    next >
Text File  |  1988-12-30  |  16KB  |  774 lines
  1. Path: xanth!nic.MR.NET!csd4.milw.wisc.edu!mailrus!ulowell!page
  2. From: page@swan.ulowell.edu (Bob Page)
  3. Newsgroups: comp.sources.amiga
  4. Subject: v02i111:  complex - complex functions
  5. Message-ID: <11028@swan.ulowell.edu>
  6. Date: 30 Dec 88 18:06:53 GMT
  7. Organization: University of Lowell, Computer Science Dept.
  8. Lines: 763
  9. Approved: page@swan.ulowell.edu
  10.  
  11. Submitted-by: kim@uts.amdahl.com (Kim E. DeVaughn)
  12. Posting-number: Volume 2, Issue 111
  13. Archive-name: applications/complex.1
  14.  
  15. Here are some complex number functions.  There are versions for both
  16. x,y and polar representations.    You know, wonderful things like +, -,
  17. *, /, sin, cos, tan, exp, ln and power.  Sorry, no hyperbolic complex
  18. functions...
  19.  
  20. Enjoy!
  21.  
  22. #    This is a shell archive.
  23. #    Remove everything above and including the cut line.
  24. #    Then run the rest of the file through sh.
  25. #----cut here-----cut here-----cut here-----cut here----#
  26. #!/bin/sh
  27. # shar:    Shell Archiver
  28. #    Run the following text with /bin/sh to create:
  29. #    complex.c
  30. # This archive created: Fri Dec 30 13:04:13 1988
  31. cat << \SHAR_EOF > complex.c
  32. /*
  33. The implementations of these functions are hereby released to the public
  34. domain.  I think it would be silly to try to sell this because:
  35. 1. limited market.
  36. 2. anyone with half a brain can do the same after looking in a frosh
  37.    level calculus textbook.
  38.  
  39. If you have find any mistakes, have a better implementation or
  40. have any questions contact me at
  41.  
  42. GEnie - R.DEVENEZIA
  43.  or
  44. 225B Hillcrest Manor Ct
  45. Utica, NY 13501
  46.  
  47.  
  48. -------------- complex math routines ---------------
  49.  
  50. Basic notation:
  51.  
  52.  z = x + iy                               one mapping       (x, y)
  53.    = r7[ cos(theta) + i sin(theta) ]      infinite mappings (r, 2pi7k7theta)
  54.    = r7cis(theta);
  55.       where cis(theta) = cos(theta) + i sin(theta)
  56.  
  57.         i7x
  58.    = r7e                                  infinite mappings (r, 2pi7k7theta)
  59.  
  60. Result Methodology:
  61.    Since most compilers can't return a structure as a function result
  62.    this implementation maintains a 'result pool'. All functions
  63.    return pointers to an entry in the result pool (set to 100 currently).
  64.  
  65.    Why do it this way? Well, I want to nest functions within functions,
  66.    and it ain't handy (hell, it ain't possible) when the result is passed
  67.    back through the argument list (you need all sorts of intermediate
  68.    variables, uggh).
  69.  
  70.    For example, I can now write an expression like
  71.           zprime = cMult (cMult (z,z), cExp (z));
  72.    voila! z' = z27exp(z)
  73.  
  74.    Dangers: Any variable that is to be a result of these functions
  75.             must be a pointer to allocated memory.
  76.  
  77.    Examples:
  78.       If predefined variables of type complex are going to be used
  79.       the addresses of these variables must be passed to the functions.
  80.  
  81.       i.e.
  82.          complex z1, z2;      / * complex variable * /
  83.          complex *result;     / * pointer to complex variable * /
  84.          real    modulus;
  85.  
  86.          z1.R = 0;
  87.          z1.i = 2;    / * 2i * /
  88.          z2.R = 1;
  89.          z2.i = 2;    / * 1 + 2i * /
  90.  
  91.          modulus = cMod (&z1);    / * argument is a pointer * /
  92.          result = cSin (&z2);     / * argument is a pointer, as is result * /
  93.          / * result is a pointer to an entry in the result pool * /
  94.  
  95.          modulus = cMod (cMult (&z1, &z2));  / * nested call * /
  96.  
  97.  
  98.       If predefined variables of type *complex are going to be used
  99.       you must first allocate memory and pass the pointer returned
  100.       [i.e. the variable] to the functions.
  101.  
  102.       i.e.
  103.          complex *z1, *z2;    / * pointers to complex variable * /
  104.          complex *result;     / * pointer to complex variable * /
  105.          real modulus;
  106.  
  107.          z1 = malloc (sizeof(complex));
  108.          z1->R = 0;
  109.          z1->i = 2;     / * 2i * /
  110.          z2->R = 1;
  111.          z2->i = 2;     / * 1 + 2i * /
  112.  
  113.          modulus = cMod (z1);    / * argument is a pointer * /
  114.          result = cSin (z2);      / * argument is a pointer, as is result * /
  115.          / * result is a pointer to an entry in the result pool * /
  116.  
  117.          modulus = cMod (cMult (z1, z2));
  118.  
  119.          and be sure to free up any memory you allocate using Amiga memory
  120.          allocation functions.
  121.  
  122.  
  123.    Personally I prefer the first method.
  124.  
  125.  
  126. available functions:
  127.  c prefix requires complex arguments
  128.  p prefix requires polar arguments
  129.  
  130.    cNextResult - maintains result pool, don't play with it
  131.    pNextResult - maintains result pool, don't play with it
  132.    cRe         - Real part
  133.    pRe         - Real part
  134.    cIm         - Imaginary part
  135.    pIm         - Imaginary part
  136.    cMod        - Modulus (radius)
  137.    pMod        - Modulus (radius)
  138.    cArg        - Argument (angle)
  139.    pArg        - Argument (angle)
  140.    cTOp        - convert complex to polar
  141.    pTOc        - convert polar to complex
  142.    cConj       - complex conjugate
  143.    pConj       - complex conjugate
  144.    MakeComplex - return a complex* from two reals
  145.    MakePolar   - return a polar* from two reals
  146.    cAdd        - addition
  147.    cSub        - subtraction
  148.    cMult       - multiplication
  149.    cDiv        - division, check for /0
  150.    cExp        - exponentiation
  151.    cLn         - natural logarithm, check for /0
  152.    cSin        - sine
  153.    cCos        - cosine
  154.    cTan        - tangent, check for /0
  155.    cIPower     - integer power
  156.    cPower      - complex to complex power
  157.    pAdd        - addition
  158.    pSub        - subtraction
  159.    pMult       - multiplication
  160.    pDiv        - division, check for /0
  161.    pExp        - exponentiation
  162.    pLn         - natural logarithm, check for /0
  163.    pSin        - sine
  164.    pCos        - cosine
  165.    pTan        - tangent, check for /0
  166.    pIPower     - integer power
  167.    pPower      - complex to complex power
  168.  
  169. */
  170.  
  171.  
  172. #ifndef sqrt
  173. #include <math.h>
  174. #endif
  175.  
  176. #define TRUE   1
  177. #define FALSE  0
  178.  
  179. #define  pi    3.1415926
  180.  
  181. #define  cPoolSize   100   /* don't put those freakin' commas */
  182. #define  pPoolSize   100   /*  after a define!! */
  183.  
  184. typedef struct {
  185.    float R;
  186.    float i;
  187. } complex;
  188.  
  189. typedef struct {
  190.    float r;
  191.    float theta;
  192. } polar;
  193.  
  194. static complex  cResultPool [cPoolSize];
  195. static polar    pResultPool [pPoolSize];
  196.  
  197. static complex  *ComplexResult;
  198. static polar    *PolarResult;
  199.  
  200. static int ComplexIndex = 0;    /* index into result pool */
  201. static int   PolarIndex = 0;    /* index into result pool */
  202.  
  203. int divisionByZero = FALSE;
  204.  
  205.  
  206. void
  207. cNextResult ()    /* arrrgh! remember, no commas here! */
  208. {
  209.    ComplexResult = &cResultPool [ComplexIndex++];
  210.    if (ComplexIndex == cPoolSize)
  211.       ComplexIndex = 0;
  212. }
  213. /* cNextResult */
  214.  
  215.  
  216. void
  217. pNextResult ()
  218. {
  219.    PolarResult = &pResultPool [PolarIndex++];
  220.    if (PolarIndex == pPoolSize)
  221.       PolarIndex = 0;
  222. }
  223. /* pNextResult */
  224.  
  225.  
  226. float
  227. cRe (z)           /* real part */
  228. complex *z;
  229. {
  230.    return (z->R);
  231. }
  232. /* cRe */
  233.  
  234.  
  235. float
  236. pRe (z)           /* real part */
  237. polar *z;
  238. {
  239.    return (z->r * cos (z->theta));
  240. }
  241. /* pRe */
  242.  
  243.  
  244. float
  245. cIm (z)           /* imaginary part */
  246. complex *z;
  247. {
  248.    return (z->i);
  249. }
  250. /* cIm */
  251.  
  252.  
  253. float
  254. pIm (z)           /* imaginary part */
  255. polar *z;
  256. {
  257.    return (z->r * sin (z->theta));
  258. }
  259. /* pIm */
  260.  
  261.  
  262. float
  263. cMod (z)          /* modulus (radius) */
  264. complex *z;
  265. {
  266.    return (sqrt ( z->R * z->R  +  z->i * z->i ));
  267. }
  268. /* cMod */
  269.  
  270.  
  271. float
  272. pMod (z)          /* modulus (radius) */
  273. polar *z;
  274. {
  275.    return (z->r);
  276. }
  277. /* pMod */
  278.  
  279.  
  280. float
  281. cArg (z)          /* angle */
  282. complex *z;
  283. /*
  284.    arg z is a one to infinite mapping. (infinite number of solutions)
  285.    arg z is a one to one mapping (one soln) using the principal determination.
  286. ie.
  287.    a + bi = (R,theta)
  288.    where theta = principal_theta 1 27pi7k; where k is an integer
  289.    and theta in (-pi,pi] is the principal determination
  290. */
  291. {
  292.    float temp;
  293.  
  294.    if (z->R == 0.0)  /* no radius */
  295.       if (z->i > 0.0) return ( pi / 2.0);
  296.       else
  297.       if (z->i < 0.0) return (-pi / 2.0);
  298.       else /* origin */
  299.          return (0.0);
  300.    else
  301.       /* Manx's atan returns the arc tangent in the range -pi/2 to pi/2 */
  302.       return (atan (z->i / z->R));
  303. }
  304. /* cArg */
  305.  
  306.  
  307. float
  308. pArg (z)          /* angle */
  309. polar *z;
  310. {
  311.    return (z->theta);
  312. }
  313. /* pArg */
  314.  
  315.  
  316. polar *
  317. cTOp (z)          /* convert complex to polar */
  318. complex *z;
  319. {
  320.    pNextResult;
  321.    PolarResult->r     = cMod (z);
  322.    PolarResult->theta = cArg (z);
  323.    return (PolarResult);
  324. }
  325. /* cTOp */
  326.  
  327.  
  328. complex *
  329. pTOc (z)          /* convert polar to complex */
  330. polar *z;
  331. {
  332.    cNextResult;
  333.    ComplexResult->R = pRe (z);
  334.    ComplexResult->i = pIm (z);
  335.    return (ComplexResult);
  336. }
  337. /* pTOc */
  338.  
  339.  
  340. complex *
  341. cConj (z)         /* complex conjugate */
  342. complex *z;
  343. {
  344.    cNextResult;
  345.    ComplexResult->R =  z->R;
  346.    ComplexResult->i = -z->i;
  347.    return (ComplexResult);
  348. }
  349. /* cConj */
  350.  
  351.  
  352. polar *
  353. pConj (z)         /* complex conjugate (polar form) */
  354. polar *z;
  355. {
  356.    pNextResult;
  357.    PolarResult->r     =    z->r;
  358.    PolarResult->theta = - (z->theta);
  359.    return (PolarResult);
  360. }
  361. /* pConj */
  362.  
  363.  
  364. complex *
  365. MakeComplex (R, i)
  366. float R, i;
  367. {
  368.    cNextResult;
  369.    ComplexResult->R = R;
  370.    ComplexResult->i = i;
  371.    return (ComplexResult);
  372. }
  373. /* MakeComplex */
  374.  
  375.  
  376. polar *
  377. MakePolar (r, theta)
  378. float r, theta;
  379. {
  380.    pNextResult;
  381.    PolarResult->r     = r;
  382.    PolarResult->theta = theta;
  383.    return (PolarResult);
  384. }
  385. /* MakePolar */
  386.  
  387.  
  388.  
  389. /*-------------------------------------------------------------------------*/
  390. /* Complex */
  391.  
  392. complex *
  393. cAdd (z1, z2)     /* addition, z1+z2 */
  394. complex *z1, *z2;
  395. {
  396.    cNextResult;
  397.    ComplexResult->R = z1->R + z2->R;
  398.    ComplexResult->i = z1->i + z2->i;
  399.    return (ComplexResult);
  400. }
  401. /* cAdd */
  402.  
  403.  
  404. complex *
  405. cSub (z1, z2)     /* subtraction, z1-z2 */
  406. complex *z1, *z2;
  407. {
  408.    cNextResult;
  409.    ComplexResult->R = z1->R - z2->R;
  410.    ComplexResult->i = z1->i - z2->i;
  411.    return (ComplexResult);
  412. }
  413. /* cSub */
  414.  
  415.  
  416. complex *         /* multiplication, z1*z2 */
  417. cMult (z1, z2)
  418. complex *z1, *z2;
  419. {
  420.    cNextResult;
  421.    ComplexResult->R = z1->R * z2->R  -  z1->i * z2->i;
  422.    ComplexResult->i = z1->R * z2->i  +  z1->i * z2->R;
  423.    return (ComplexResult);
  424. }
  425. /* cMult */
  426.  
  427.  
  428. complex *
  429. cDiv (z1, z2)     /* division, z1/z2 */
  430. complex *z1, *z2;
  431. {
  432.    float z2_modulusSquared;
  433.  
  434.    cNextResult;
  435.    z2_modulusSquared = z2->R * z2->R  +  z2->i * z2->i;
  436.    if (z2_modulusSquared == 0.0) { 
  437.       divisionByZero = TRUE;
  438.       ComplexResult->R = 0.0;
  439.       ComplexResult->i = 0.0;
  440.    }
  441.    else {
  442.       divisionByZero = FALSE;
  443.       ComplexResult->R = (z1->i * z2->i  +  z1->R * z2->R)
  444.                 / z2_modulusSquared;
  445.       ComplexResult->i = (z1->i * z2->R  -  z1->R * z2->i)
  446.                 / z2_modulusSquared;
  447.    }
  448.    return (ComplexResult);
  449. }
  450. /* cDiv */
  451.  
  452.  
  453. complex *
  454. cExp (z)          /* exponentiation */
  455. complex *z;
  456. /*
  457.           infinity
  458.     z        --   n
  459.    e  == 1 + \   z
  460.              /  ----
  461.              --  n!
  462.             n=1
  463.  
  464. OR the method used (by the following equality):
  465.  
  466.     z     x + iy     x   iy      x                x
  467.    e  == e       == e 7 e   ==  e * cos(y) + i * e * sin(y) 
  468.  
  469. */
  470. {
  471.    float eToTheX;
  472.  
  473.    cNextResult;
  474.    eToTheX = exp (z->R);
  475.    ComplexResult->R = eToTheX * cos (z->i);
  476.    ComplexResult->i = eToTheX * sin (z->i);
  477.    return (ComplexResult);
  478. }
  479. /* cExp */
  480.  
  481.  
  482. complex *
  483. cLn (z)           /* natural logarithm */
  484. complex *z;
  485. /*
  486.    ln z is a one to infinite mapping. (infinite number of solutions)
  487.    ln z is a one to one mapping (one soln) using the principal determination.
  488.  
  489.    complex   real
  490.        ln z =  ln(|z|) + i * arg(z)
  491. */
  492. {
  493.    cNextResult;
  494.    if (cMod (z) == 0.0) {
  495.       divisionByZero = TRUE;
  496.       ComplexResult->R = 0.0;
  497.       ComplexResult->i = 0.0;
  498.    }
  499.    else {
  500.       divisionByZero = FALSE;
  501.       ComplexResult->R = log ( cMod (z) );
  502.       ComplexResult->i = cArg (z);   /* principal determination */
  503.    }
  504.  
  505.    return (ComplexResult);
  506. }
  507. /* cLn */
  508.  
  509.  
  510. complex *
  511. cSin (z)          /* sine */
  512. complex *z;
  513. /*
  514.           infinity
  515.              --        k
  516.    sin z  == \     (-1)     2k+1
  517.              /  ---------- z
  518.              --  (2k + 1)!
  519.              k=0
  520.  
  521. OR the method used (by the following equality):
  522.  
  523.              1     iz    -iz
  524.    sin z == --- [ e   - e   ]
  525.             27i
  526. also,
  527.          == sin(x)7cosh(y) + i cos(x)7sinh(y)
  528.               where cosh is the hyperbolic cosine and
  529.                     sinh is the hyperbolic sine.
  530. */
  531. {
  532. return (cMult (MakeComplex (0.0, -0.5),   /* 1/27i */
  533.                cSub (cExp (cMult (MakeComplex (0.0, 1.0), z)),
  534.                      cExp (cMult (MakeComplex (0.0,-1.0), z))
  535.                     )
  536.               )
  537.        );
  538. }
  539. /* cSin */
  540.  
  541.  
  542. complex *
  543. cCos (z)          /* cosine */
  544. complex *z;
  545. /*
  546.           infinity
  547.              --      k
  548.    cos z  == \   (-1)   2k
  549.              /  ------ z
  550.              --  (2k)!
  551.              k=0
  552.  
  553. OR the method used (by the following equality):
  554.  
  555.              1     iz    -iz
  556.    cos z == --- [ e   + e   ]
  557.              2
  558. also,
  559.          == cos(x)7cosh(y) - i sin(x)7sinh(y)
  560.               where cosh is the hyperbolic cosine and
  561.                     sinh is the hyperbolic sine.
  562. */
  563. {
  564. return (cMult (MakeComplex (0.5, 0.0),
  565.                cAdd (cExp (cMult (MakeComplex (0.0, 1.0), z)),
  566.                      cExp (cMult (MakeComplex (0.0,-1.0), z))
  567.                     )
  568.               )
  569.        );
  570. }
  571. /* cCos */
  572.  
  573.  
  574. complex *
  575. cTan (z)          /* tangent */
  576. complex *z;
  577. /*
  578.              sin z
  579.     tan z = -------
  580.              cos z
  581. */
  582. {
  583.    if (cMod (z) == 0.0) { 
  584.       divisionByZero = TRUE;
  585.       cNextResult;
  586.       ComplexResult->R = 0.0;
  587.       ComplexResult->i = 0.0;
  588.       return (ComplexResult);
  589.    }
  590.    else {
  591.       divisionByZero = FALSE;
  592.       return (cDiv (cSin (z), cCos (z)));
  593.    }
  594. }
  595. /* cTan */
  596.  
  597.  
  598. complex *
  599. cIPower (z, n)    /* complex number to an integer power */
  600. complex *z;
  601. int      n;
  602. /*
  603.     n      n
  604.    z  = |z| 7 cis (n7theta)
  605. */
  606. {
  607.    float modulusToTheN, ntheta;
  608.  
  609.    cNextResult;
  610.    modulusToTheN = pow (cMod (z), n);        /* Manx's power function */
  611.    ntheta = n * cArg (z);
  612.    ComplexResult->R = modulusToTheN * cos (ntheta);
  613.    ComplexResult->i = modulusToTheN * sin (ntheta);
  614.    return (ComplexResult);
  615. }
  616. /* cIPower */
  617.  
  618.  
  619. complex *
  620. cPower (z1, z2)   /* complex number raised to a complex number */
  621. complex *z1, *z2;
  622. /*
  623.      z2      z2*log(z1)
  624.    z1   = exp
  625.                                                                        z2
  626. Note: There are an infinite number of points, (r,theta), which equal z1  .
  627.       There is only one point if we take the principal determinations of
  628.       log and exp.
  629. */
  630. {
  631.    return (cExp (cMult (z2, cLn (z1))));
  632. }
  633. /* cPower */
  634.  
  635.  
  636. /*-------------------------------------------------------------------------*/
  637. /* Polar */
  638.  
  639. polar *
  640. pAdd (z1, z2)     /* addition */
  641. polar *z1, *z2;
  642. {
  643.    return (cTOp (MakeComplex (pRe (z1) + pRe (z2),
  644.                               pIm (z1) + pIm (z2))
  645.                 )
  646.           );
  647. }
  648. /* pAdd */
  649.  
  650.  
  651. polar *
  652. pSub (z1, z2)     /* subtraction */
  653. polar *z1, *z2;
  654. {
  655.    return (cTOp (MakeComplex (pRe (z1) - pRe (z2),
  656.                               pIm (z1) - pIm (z2))
  657.                 )
  658.           );
  659. }
  660. /* pSub */
  661.  
  662.  
  663. polar *
  664. pMult (z1, z2)    /* multiplication */
  665. polar *z1, *z2;
  666. {
  667.    float theta;
  668.  
  669.    pNextResult;
  670.    PolarResult->r     = z1->r * z2->r;
  671.    theta = z1->theta + z2->theta;
  672. /*
  673.    principal determination ?
  674.    if ((theta <= -pi) OR (theta > pi)) 
  675.       theta = theta + 2*pi* INT(theta/2/pi)
  676.    }
  677. */
  678.    PolarResult->theta = theta;
  679.    return (PolarResult);
  680. }
  681. /* pMult */
  682.  
  683.  
  684. polar *
  685. pDiv (z1, z2)     /* division */
  686. polar *z1, *z2;
  687. {
  688.    pNextResult;
  689.    if (z2->r != 0.0) {
  690.       divisionByZero = FALSE;
  691.       PolarResult->r = z1->r / z2->r;
  692.    }
  693.    else {
  694.       divisionByZero = TRUE;
  695.       PolarResult->r = 0.0;
  696.    }
  697.    PolarResult->theta = z1->theta - z2->theta;
  698.    return (PolarResult);
  699. }
  700. /* pDiv */
  701.  
  702.  
  703. polar *
  704. pExp (z)          /* exponentiation */
  705. polar *z;
  706. {
  707.    return (cTOp (cExp (pTOc (z))));
  708. }
  709. /* pExp */
  710.  
  711.  
  712. polar *
  713. pLn (z)           /* natural logarithm */
  714. polar *z;
  715. {
  716.    return (cTOp (cLn (pTOc (z))));
  717. }
  718. /* pLn */
  719.  
  720.  
  721. polar *
  722. pSin (z)          /* sine */
  723. polar *z;
  724. {
  725.    return (cTOp (cSin (pTOc (z))));
  726. }
  727. /* pSin */
  728.  
  729.  
  730. polar *
  731. pCos (z)          /* cosine */
  732. polar *z;
  733. {
  734.    return (cTOp (cCos (pTOc (z))));
  735. }
  736. /* pCos */
  737.  
  738.  
  739. polar *
  740. pTan (z)          /* tangent */
  741. polar *z;
  742. {
  743.    return (cTOp (cTan (pTOc (z))));
  744. }
  745. /* pTan */
  746.  
  747.  
  748. polar *
  749. pIPower (z, n)    /* polar number to an integer power */
  750. polar *z;
  751. int    n;
  752. {
  753.    pNextResult;
  754.    PolarResult->r     = n * z->r;
  755.    PolarResult->theta = n * z->theta;
  756.    return (PolarResult);
  757. }
  758. /* pIPower */
  759.  
  760.  
  761. polar *
  762. pPower (z1, z2)   /* polar number raised to a polar number */
  763. polar *z1, *z2;
  764. {
  765.    return (cTOp (cPower (pTOc (z1), pTOc (z2))));
  766. }
  767. /* pPower */
  768. SHAR_EOF
  769. #    End of shell archive
  770. exit 0
  771. -- 
  772. Bob Page, U of Lowell CS Dept.  page@swan.ulowell.edu  ulowell!page
  773. Have five nice days.
  774.